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ABSTRACT 


Performance models of the Space Shuttle Main Engine (SSME) 
contain iterative strategies for determining approximate 
solutions to nonlinear equations reflecting fundamental mass, 
energy, and pressure balances within engine flow systems. Both 
univariate and multivariate Newton-Raphson algorithms are 
employed in the current version of the engine Test Information 
Program (TIP). Computational efficiency and reliability of these 
procedures is examined. A modified trust region form of the 
multivariate Newton-Raphson method is implemented and shown to be 
superior for off nominal engine performance predictions. A 
heuristic form of Broyden's Rank One method is also tested and 
favorable results based on this algorithm are presented. 
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INTRODUCTION 


Predictions of the steady-state operational characteristics 
of the Rocketdyne Space Shuttle Main Engine (SSME) are provided 
by computer programs which model and analyze engine sytem 
performance. The Test Information Program (TIP) is a FORTRAN 
based engine analysis software package which performs the 
following functions. 

1. Power Balance - Simulates engine performance by balancing 

mass and energy flows for assumed nominal 
operation of SSME components. 

2. Data Reduction - Uses actual test data to define operating 

characteristics of a specific SSME by 
adjusting component performance parameters. 

3. Base Balance - Refines operating predictions of the Data 

Reduction portion by adjusting nine 
performance variables. 

4. Rated Portion - Extends refined performance model to 

simulate actual engine operation over a 
range of conditions. 

TIP balances mass and energy flows during engine performance 
prediction, and balances theoretical results with test 
information to refine these predictions. Figure 1 displays the 
SSME flow system network that is balanced by TIP to ensure 
satisfaction of the conservation of mass and conservation of 
energy principles as well as adherence to limitations imposed by 
the Second Law of Thermodynamics. 

Since flow processes within the SSME are governed by a set 
of nonlinear equations, iterative techniques are required to 
computationally predict a balanced steady-state flow condition. 
Two subroutines within the TIP code perform iterative nonlinear 
equation solving functions. These routines are described below. 

TIP Iteration Subroutines 

1. NLEST - a single nonlinear equation root finding routine 
CALL statement accessed as subroutine NLEST 
ENTRY statement accessed as subroutine NLREST 

NLEST - Find z such that (x-y) = F(z) = 0 . 

NLREST - Find x such that (x-y) = F(x) = 0 
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2. NLESTN - a simultaneous nonlinear equation solver 

CALL statement accessed as subroutine NLESTN 

NLESTN - Find x = (xl,x2 . . . xN) such that 

yl = Fl(x) = 0 
y2 = F2( x) = 0 


yN = FN( x) = 0 

The basic function of TIP is to provide steady-state 
operation and performance predictions. Engine flow conditions 
are established by solving systems of nonlinear balance equations 
arising from fundamental flow mechanics. Hence, use of the above 
described solver routines is fundamental to and pervasive within 
TIP aj evidenced by the ninty-four calls to NLEST (or NLREST) and 
NLESTN which occur in the current version of this performance 
analysis package. 

In addition to being pervasive, the iteration loops 
initiated by calls to subroutines NLEST and NLESTN are programmed 
in a complex sequence with considerable nesting and crossover. 

As an example of this complexity, consider the iteration loop 
sequencing of subroutine BAL exhibited in Appendix 1. Subroutine 
BAL is the routing routine for the Power Balance component of 
TIP. Significant iteration loop nesting, as many as seven loops 
deep, and loop crossover are apparent upon examination of this 
listing . 


OBJECTIVES 


Because of the intrinsic importance of the iterative 
nonlinear equation solvers in the TIP code, it is important that 
these procedures be both reliable and efficient. The objectives 
of this research effort are 

1. To evaluate the iterative schemes employed within the TIP 
performance model 

2. To modify and test these schemes as suitable to achieve greater 
reliability and efficiency 

3. To perform a cursory review of fundamental TIP code logic and 
procedure 
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NONLINEAR EQUATION SOLVERS 


The nonlinear equation root finding scheme currently 
employed in the TIP subroutine NLEST is a heuristic procedure 
which combines a secant method searoh with a false positions 
method (see e.g. Burden and Faires [1]). Initially the secant 
method is used to update the root approximation. Whenever 
successive applications of the secant method establish a bracket 
about the root, NLEST switches to a false positions technique. 

It is possible to lose the established root bracket due to loop 
nesting and the accompanying effects of shifting due to inner 
loop convergence tolerances. Because of this an ad hoc strategy 
for switching between secant and false positions searches is 
employed . 

Within the context of the TIP code, the existing root 
finding scheme in NLEST was tested against several other well 
known procedures. The procedures implemented and tested by this 
author are listed below. 

1. Finite difference Newton-Raphson 

2. Secant (only) 

3. False positions with step bracketing 

4. Quadratic interpolation polynomial 

In few instances were results superior to those obtained using 
the current NLEST algorithm, and in many cases the overall 
performance and convergence of the TIP code were adversely 
affected. Hence, the heuristic root finding strategy implemented 
within subroutine NLEST appears to be a robust and effective 
method within the TIP performance analysis model. 

The simultaneous nonlinear equation solver currently 
implemented within subroutine NLESTN is a finite difference 
multivariate Newton-Raphson method (often referred to as simply 
Newton's method, see e.g. reference [1], page 496) with a 
direction skewing trust region boundary. To evaluate the 
effectiveness of this algorithm it is necessary to refer to the 
theoretical basis of the Newton-Raphson method. This basis will 
be outlined below. 



The fundamental problem addressed by simultaneous nonlinear 
equation solvers can be expressed mathematically as follows. 

Find x = ( xl, x2 . . . xN ) T in R* 

such that F = ( Fl(x), F2(x) . . . FN(x) ) T = 0 


or equivalently Fl(x) = 0 

F2( x) = 0 


( 1 ) 


where 


FN ( x) = 0 

F : R N 


At stage k the fundamental iteration of the Newton-Raphson 
method for solving problems of this type can be written 

x k+l = x k ~ Jk 1 * Fk (2) 

where 

Ffc = F(x k ) 

and 

J k - the Jacobian of F evaluated at x k 


Rapid convergence of the Newton-Raphson method to the solution of 
(1) can be established rigorously if, for some k, x k is 

sufficiently close to the solution vector designated x*. 

Formal conditions for convergence are stated in the following 
theorem (see e.g. Rheinboldt [2] ). 

Theorem . 

If F has continuous first partial derivatives in some 
neighborhood of t he solution x , if the Jacobian of F is 
nonsingular at x and its elements satisfy a Lipschitz condition, 
and if x k is sufficiently close to x for some k, then the 
Newton- Raphson method is well defined for all k and converges at 
second order, i.e. there exits a positive integer m and a 
positive real number b such that 

||f k+ 1 !| / ||F k |J 2 < b whenever k > m 


Despite the second order (or quadratic) convergence 
indicated in the above theorem, the Newton-Raphson method suffers 
from two serious disadvantages from the point of view of 
practical calculation. First, computation of the Jacobian matrix 
at each stage of the iteration is extremely costly in terms of 
computer resources. Often analytical partial derivatives are not 
available and finite difference approximations lack the precision 
necessary for ultimate convergence. The second disadvantage of 
the Newton-Raphson method arises from the need to have a 
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sufficiently accurate initial estimate of the solution in order 
to guarantee convergence. Satisfaction of this requirement is 
impossible to measure before initiation of the iteration sequence 
and often difficult to obtain in practice. 

Extension of the basic Newton-Raphson procedure to include a 
subiteration or line search have been somewhat successful in 
removing the accurate initial estimate requirement. These 
methods include a strategy to select a positive real number 
lambda such that the iteration scheme 

*k+l = *k + lambda * J^ 1 * F k (3) 

reduces some measure of error, typically 


ll'k+lll < l|Fk|| 


<4 > 


at each stage. Although Newton-Raphson methods incorporating 
line searches extend the domain of convergence, they do so with 
significant computational overhead. 

Trust region or restricted step methods are a compromise 
between convergence limited self-scaling iteration procedures and 
computationally intense methods incorporating line searches. 

These methods simply provide an upper bound on the distance 
between iteration steps. This bound may be absolute or scaled by 
position within the domain of x. 

The current version of subroutine NLESTN provides a 
multivariate Newton-Raphson method with a trust region approach. 
Unfortunately, the trust region bound is applied componentwise on 
x which has the effect of skewing the correct Newton-Raphson 
method search direction. This skewing process removes any 
theoretical convergence characteristics and indicates the 
possibility of convergence difficulties when the iteration 
procedure is initiated at a point remote from the immediate 
vicinity of the solution. Difficulties of this type are 
currently experienced as will be discussed in the next section. 

In an effort to correct the skewing problem inherent in the 
current version of NLESTN, a modified trust region form of the 
multivariate Newton-Raphson method was implemented within TIP. 

The basic iteration sequence of this method is given by 


where 


x k+l = *k 


lambda = min 


- lambda * * F k 

|| PCTMAX * x k || 
l|jk 1 *F k j| 


(5) 

( 6 ) 


PCTMAX is a user defined parameter which serves to scale the size 
of the trust region. A comparison of search steps obtained 
using the skewed direction trust region approach and the 
corrected strategy are displayed in Figure 2. A complete listing 
of the modified trust region form of NLESTN implemented within 
TIP is given in Appendix 2. Results using this implementation 
are presented in the next section. 

A trust region form of the Broyden Rank One [3] nonlinear 
system solver was also implemented and tested in an effort to 
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Skewed x k+l 


Figure 2. Comparison of skewed trust region search vector and 
correct trust region search vector. 
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reduce the computational overhead associated with Jacobian 
calculations and linear system solvers in the Newton-Raphson 
algorithm . The Broyden Rank One method is a modified Newton- 
Raphson iteration in which the Jacobian matrix is replaced by an 
approximation B which is updated at each iteration by using 
current information about F. In the Broyden strategy, a rank one 
modification of B is made at each iteration. The basic method is 
described below. 

Basic Broyden Rank One Method 

0. Select x Q the iteration starting point and 

B 0 an N x N matrix approximation to the inverse 
Jacobian of F at x D . 

For k = 0,1,2... 

1. Evaluate 

2 . Set s k = “ B k * F k 

3. Set x k+1 = x k + s k 

4. Evaluate *k + l 

5. Set y k = F k+1 - F k 

6. If y k = 0 Set B k+ i = B k 

Otherwise set B k+1 = B k + <s k - B k y k ) * v T 

where v T = s T * B k / (s T * B k * yjj) 

Theoretical convergence properties of the basic Broyden Rank 
One algorithm have been clarified by Gay [4,5]. Although not as 
rapidly convergent theoretically as the Newton-Raphson method, 
testing has often shown the Broyden technique to be superior [3] 

A complete listing of the Broyden Rank One trust region form 
of NLESTN implemented within TIP is given in Appendix 3. Results 
using this implementation are presented in the next section. 

Function minimization strategies have often been used to 
solve nonlinear systems. A class of specialized minimization 
techniques (see e.g. [6]) been developed which address least 
squares problems of the type 

minimize F(x) T * F(x) by selection of x (7) 

General nonlinear minimization algorithms can also be employed to 
solve problems such as (7) which arise naturally from nonlinear 
equation systems. These methods are typically computer intensive 
although potentially effective and robust. 

Although no testing of TIP performance using sophisticated 
minimization methods was conducted, such techniques hold promise 
lor future application within engine performance models. 
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COMPUTATIONAL RESULTS 


Computations involving the TIP performance model were 
conducted on the EADS network at Marshall Space Flight Center. 
An IBM 3084 was the host processor. Within TIP, only Power 
Balance model prediction analyses were run. Initializing input 
data was fixed for all computations with the exception of 
LOX-fuel ratio and power level expressed as a percentage of 
engine rated power. Values for these flow characteristics were 
specified within and outside of the nominal operating ranges 
defined below. 


Nominal SSME Operating Range 

LOX-Fuel ratio 6.0 - 6.5 

Power Level (%) 65 - 109 

Variation of these operating control parameters permitted an 
evaluation of the function of the nonlinear equation solving 
algorithms over a range of practical problems. 

Output from a typical Power Balance model analysis is 
extensive and includes many operating and performance 
characteristics of the SSME subsystems. These physical data will 
not be presented. Since the object of this investigation was to 
test and evaluate performance model iteration methods, 
computational performance data involving subroutines NLEST and 
NLESTN were collected. These data were organized and presented 
for each convergent Power Balance analysis as shown in Table 1. 
The prescribed values of mixture ratio and power level for the 
specific analysis are displayed at the top of this table. In 
addition, the univariate and multivariate iteration procedures 
used in the particular analysis are presented. Performance 
evaluation data for each iterative procedure are presented as 
described below. 


NLEST Loop Summary Information 


Loop 

Entered 

Closed 

Max 

Average 


iteration loop identification number 
number of times specific iteration loop was entered 
number of times specific iteration loop successfully 
converged and closed 

maximum number of iterations required for convergence 
of specific iteration loop 

average number of iterations required for convergence 
of specific iteration loop 
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Table 1. Iteration loop performance data for a typical Power 
Balance model analysis. 


MIXTURE RATIO - 6.5 

NLEST - ORIGINAL 
NLESTN - BROYDEN RANK ONE 

NLEST LOOP SUMMARY LOOP 

1 

2 

3 

5 

6 

7 

8 
9 

IX 

12 

14 

15 

16 

17 

18 

19 

20 
21 

23 

24 

25 

26 

27 

28 

30 

31 

32 

35 

36 

3 7 

38 

39 

40 

42 

43 

4 5 
46 

48 

49 

50 
60 
61 
62 
67 
70 

NLESTN SUMMARY LOOP ENTER 

1 510 

2 19 

3 769 

4 294 

5 286 

6 3375 


POWER LEVEL - 109 % 


ENTERED 

71 

67 

589 

6 

10 

434 

315 

87 

149 

788 

1002 

550 

1029 

119 

1100 

657 

679 

13 

1113 

2108 

861 

1698 

1737 

977 

117 

5528 

191 

442 

3 

5 

429 

191 

5 

856 

69 

69 

769 

388 

431 

174 

594 

177 

122 

3 

4 


CLOSED 

67 

59 

294 

3 

5 

428 

174 

69 

67 

428 

550 

294 

510 

69 

428 

290 

290 

5 

428 

720 

428 

788 

856 

428 

69 

2517 

69 

69 

1 

1 

428 

191 

1 

428 

19 

69 

408 

294 

294 

87 

360 

122 

67 

1 

1 


JACOB CLOSE MAX 
157 149 5 

a 1 8 

361 136 3 

6 5 69 3 

79 69 3 

675 675 1 



JE/C 
1.054 
8.000 
2.654 
0 .942 
1.145 
1.000 


AVERAGE 
1 . 060 
1.136 
2.003 
2.000 
2.000 
1.014 
1 . 810 
1.261 
2.224 
1 .841 
1.822 
1 . 871 
2.018 
1.725 
2 .570 
2 . 266 
2 . 341 
2.600 
2 . 600 
2.928 
2.012 
2 .155 
2 . 029 
2.263 
1.696 
2.196 
2 . 768 
6.406 

3 . 000 
5 .000 
1.002 
1 . 000 

5.000 
2 .000 
3.632 
1 .000 
1.885 
1.320 
1.466 

2.000 
1.650 
1.451 
1 .821 

3 .000 

4.000 

LP/C 
3 .423 
19.000 
5.654 
4 .261 
4.145 
5.000 
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NLESTN Loop Summary Information 


Loop 

Enter 

Jacob 

Close 

Max 

JE/C 

LP/C 


same as above 

same as above for entered 

total number of Jacobian or Jacobian approximation 
evaluations 

same as above for closed 
same as above 

average number of Jacobian evaluations per convergent 
iteration 

average number of iteration loop passes per 
convergent iteration 


The amount of effort expended in iteration processes is 
evident from Table 1, with over 15,000 univariate loop entries 
and over 5,000 multivariate loop entries documented for this 
specific Power Balance analysis in order to achieve convergence. 
This effort level is typical of convergent analyses performed in 
this study. 

In order to compare the efficiency and reliability of 
various NLESTN implemented multivariate nonlinear equation 
solving methods, the total number of NLESTN loop entries was 
tabulated for each of several analyses using different 
multivariate iteration strategies. These data are presented in 
Table 2. Examination of the information in Table 2 suggests the 
following 

Conclusions Based on Table 2 Data 


1. The current (N-R Orig) Newton-Raphson implementation often 
fails to converge for mixture ratios or power levels outside 
the nominal region. This result was expected since the 
direction skewing trust region method forces the iteration 
sequence to take less appropriate steps for conditions that 
cause the solution to be further removed from the initiation 
data. It is notable that within the nominal region for 
mixture ratio and power level, the current Newton-Raphson 
implementation performs almost the same as the modified 
method (N-R Mod). This occurs because within the nominal 
operating range, the trust region boundary is not reached 
during the search process since the initiation data is close 
to the converged solution. 

2. The modified Newton-Raphson method (N-R Mod) with corrected 
trust region is more efficient, requiring fewer total NLESTN 
loop passes, for conditions outside and on the boundary of 
the nominal regions for mixture ratio and power level. This 
is due to the corrected search direction method employed at 
the trust region boundary. In addition, the modified Newton- 
Raphson method is more reliable than the original method, 
converging for several cases with outside nominal mixture 
ratios or power levels. 
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3 . The Broyden Rank One method implementation is extremely 
ia 5J e ’ , having converged for all cages considered. 

The efficiency of this method is somewhat erratic, often 
requiring substantially fewer loop passes for convergence than 
he Newton Raphson methods and yet occasionally requiring 

* t “^?“ y h ;S re . 9ffort t0 « * converged^solution . 

This erratic behavior was not wholly unexpected due to the 

P H r ? X i 10 u te ? ature ^he Jacobian estimate employed and 
updated by the algorithm. H y a 

Further comparisons of the modified Newton-Raphson and 

a represented^ in Tables "s' 1 ‘f! ourront SLEST * i«Ple«.entation 
are presented m Table 3. Results are presented only for 

nalyses m which the current multivariate iteration strategy 

the^mod if ied V Newr nC6 R I h * in,pr ° ved efficiency gained by use of 
the modified Newton-Raphson method outside the nominal operating 

- a *? ln evident. The erratic efficiency of the B?oyden 
method is clearly displayed. y 
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Table 2. 


Loop passes through multivariate subroutine NLESTN. 


MR 

PL(%) 

M-R Orig 

N-R Mod 

Broyden 

o 

CD 

65 

9,077 

9,077 

2,527 

6.0 

100 

10,586 

10,586 

6,099 

6.011 

104 

4,993 

4,993 

14,711 

6.5 

104 

9,152 

9,152 

6,515 

6.5 

109 

13,308 

7,209 

5,253 

6.6 

109 

F 

11,491 

8,009 

6.7 

109 

F 

11,514 

13,434* 

6.8 

109 

F 

24,935* 

10,303* 

6.5 

112 

11,757 

11,723 

13,309 

6.5 

115 

15,118 

12,213 

21,829 

6.5 

120 

F 

F 

7,809 


F - Failure to converge to specified tolerance in allowed number 
of iterations 

* - Trust region interval reduced to +- 5% of current independent 
variable value 


Table 3. Change (%) in number of loop passes through subroutine 
NLESTN using the original Newton-Raphson method results 
as standard. 


MR 

PL(X) 

6.0 

65 

6 . 0 

100 

6.011 

104 

6.5 

104 

6.5 

109 

6 . 5 

112 

6 . 5 

115 
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% Change 


R Mod 

Broyden 

0 

- 72.2 

0 

- 42.4 

0 

+194.6 

0 

- 28.8 

45.8 

- 60.5 

0.3 

+ 13.2 


19.2 


+ 44.4 


RECOMMEND AT I OHS 


The following recommendations for improvement of the 
iterative procedures within the Test Information Program are 
motivated by the computational results described in the previous 
section of this report and by investigation of TIP logic. 

Recommendations for Improvement of Iterative Procedures 

1. Immediately implement modified multivariate Newton-Raphson 
method with corrected trust region approach in subroutine 
NLESTN . The subroutine described in Appendix 2 is one 
implementation of this method. 

2 . Continue to test and refine the potentially effective Broyden 
Rank One method for the iterative solution of simultaneous 
nonlinear equations. The subroutine described in Appendix 3 
is one implementation of this method. 

3. Perform computational experimentation using flexible loop 
tolerances and flexible trust region bounds in the iterative 
routines. These modifications could substantially improve 
the efficiency of the TIP iteration sequence. 

4. Incorporate and test a formal line search algorithm within the 
multivariate iteration scheme to enhance convergence and 
reliability for strongly off nominal engine operation. 

5. Perform a detailed sensitivity analysis of all iteration loop 
independent variables to determine uncertainty limits 
associated with loop tolerances. 

6. Review iteration loop logic sequencing. Modify and test 
sequencing to achieve improved computational efficiency. 

In addition to the recommendations involving computational 
procedure listed above, a limited review of the TIP code 
motivates the structuring recommendations provided below. 

Code Structuring Recommendations 

1 Clearly identify and separate TIP program components. 
Theoretical base (flow physics) 

Computational base (formal numerical algorithms) 
Experimental base (engineering performance parameters and 

other approximations) 
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2. Clearly identify the following. 

Independent variables of model analysis (user defined and 
controlled physical inputs) 

Arbitrarily prescribed and constant parameters of model 

analysis (code designer defined and restricted input) 
Dependent variables of model analysis (solution variables 
requiring initial approximation) 

3. Review prescribed flow and performance variable dependencies 
in model for accuracy and completeness. 

4 ‘ an organized data input structure descriptive of 

SSME flow systems, i.e. number nodes, branches, and devices 
and formally identify connectivity within the data structure. 

3. Fully document TIP program physical logic sequence. 

6. Construct postprocessors that clearly exhibit physical 
balances for appropriate engine subsystems as a means of 

verification . 


These recommendations 
improve confidence in 


are very basic and if implemented will 
and reliability of TIP analysis results. 
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Appendix 2. Subroutine NLESTN implementation with corrected 
trust region Newton-Raphson method. 


*********** 


********** 


CODE DESIGNATION TRUST 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


r 

c * 


c 


3 KUUTINE NLESTN <ID,N,A,KK,X1,Y1,T1,F1 X2 Y2 T2 F 2 vi v * 

+ X4 , Y4 r T4 ,F4,X5,Y5,T5,F5,X6,Y6 f T6 r F6 ) ' ' T2 ' F2 ' X3 ' Y3 ^3,F3 ( 

IMPLICIT REAL * 8 ( A-H , 0 - Z ) 


CALL STATEMENT . . . 

CALL NLESTN { I D „ N , A , KK , Xl , Y 1 , T1 , F 1 , X 2 , Y 2 T2 F2 
-••WHERE N IS ).6. ' ' 2 * ‘ 


. XN , YN , TN , FN ) 


elpi-ssssr 

« ;;'rL"™o» s ”'" E ” s,oH '° « ™ '«««««. 

SINGULAR MATRIX. 

ALL^S^Rr^fcc ° F ITERATI0NS “AS BEEN EXCEEDED. 

I„ S ARE LESS THAN tolerance, iteration complete 

ONE OR MORE Y * S ARE GREATER THAN TOLERANCE , RE ITERATE . 


KK — - 3 
KK = - I 
KK= 0 
KK= I 


INCLUDE (INSAVE) 
INCLUDE (DPSAVE) 


DIMENSION B(36) , SX ( 6 ) , S Y ( 6 ) , ST ( 6 ) 
DIMENS ION D;6 # 6) , A ( N , N ) ,C(36> 

DATA LIMIT , PCTMAX/20 , 0 . 2/ 


, SF ( 6 ) , SDF ( 6 ) 


sx« 1 .! = l 
s x ( 2 ) = x 2 
s:<: i 3 * - X j 

SX f 5 )*v5 
SX < 6 } =s x r, 


ARRAY INITIALIZE SUBROUTINE INPUT ARGUMENTS 


SY ( 1 ) = Y 1 
SY ( 2 ) - Y 2 
SY( 3 ) = Y 3 
SY( 4 ) = Y 4 
SY( 5)=Y5 
SY( 6 )=Y 6 


DRfCiNA? 

Rooiv 


WALfiy 
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ORIGINAL PAGE S3 
OF POOR QUALITY 


c 


c 

c * * 


10 


STM i )-Tl 
ST{ 2 )=T2 
ST{ 3)=T3 
ST ( 4 )»T4 
ST( 5 ) »T5 
ST ( 6 )=T 6 

SF(I)*F1 
SF ( 2 }*F2 
SF( 3 )=F3 
5 F ( 4 ) »F 4 
SF(5)=F5 
SF ( 6 ) = F 6 


L=LOOPN( ID ) 

NUMN ( ID f 1 ) aNUMN ( ID, 1 ) +1 
I CONV = 0 
DO 10 1 = 1, N 

IF (ICONV.GT.O) GO TO 10 

IF ( DABS ( S Y ( I ) ) . GT . DABS ( ST { I ) ) ) ICONV-1 

CONTINUE 

IF UCONV.GT.O) GO TO 20 
NUMN ( ID, 3 )=NUMN< ID, 3 }+l 

NUMN ( I D , 4 ) =MAX 0 ( NUMN ( ID , 4 ) , KOUNL ( ID ) ) 

KOUNL ( ID) =0 

KK=0 

GO TO 200 


STEP 1 , 


CHECK FOR CONVERGENCE 


C 

C * * 
20 

901 


FORMAT ( J ' 901 ] LIMIT ' ID » ( S Y( 1 ) ,ST< I ) , SX< I ) , 1=1 ,N) 


+ 1H 
+ 1H 
+ 1H 
+ 

K = i 


.'ERROR IN NLESTN , NO SOLUTION WITHIN I 3 ITERATIONS ■ 
FIRST ARGUMENT IN THE CALL STATEMENT IS ',15 / 

'ERROR VALUES TOLERANCES INDEPENDENT VARIABLES' / 

3 ( ill 5 . 6 ) ) " ’ 




GO TO 130 



30 DO 40 I =1 , n 
40 SDF( I ) =SF ( I ) — 1 . 0 

IF ( L . GT . 0 ) GO TO 6 0 
L=0 

DO 50 1 = 1, n 
A(I,N)=SY(I) 

50 A { I , N- 1 ) =SX { I ) 

KO ON L ( I D ) =KOUNL ( ID) +1 
NUMN( ID , 2 ) =NUMN { ID, 2 ) +1 
GO TO 180 


oTEP 3. INITIALIZE ITERATION SEQUENCE 


60 


IF ( L . LT . N ) GO TO 8 0 


STEP 4. 


ESTIMATE JACOB’IAN PARTIALS 
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70 

80 


100 

c 

c * * * 

c 

c 

c 

c 

c 


110 


c 

c * * i 


902 


DO 70 I = 1 , N 
B ( I ) =A( I , N } 

SX ( L )*SX< L )/SF { L ) 
DP*SX(L) * S DF ( L ) 

DO 100 1 = 1, N 

A(I,L)=(SY(I) —A ( I , N ) ) /DP 
IF (L.LT.N) GO TO 180 


ORfGiM/iL paw 
OF PCCfc O 


;tv 


****************** 


STEP 5. 


DO 110 1=1 

DO 110 J=1 
D( I , J)=A( I 
SCALE-0 . 
K-ISIMDD ( 6 , N 


N 

N 

J) 


DETERMINE NEWTON-RAPHSON STEP INCREMENT 
COLUMN 1 OF ARRAY D AS RETURNED FROM 
FUNCTION ROUTINE ISIMDD IS THE DOT 
PRODUCT OF THE INVERSE JACOBIAN WITH 
THE DEPENDENT VARIABLE VECTOR 
( Y1,Y2. . .YN) 


1 , D, B, SCALE , C) 


IF { K. 

WRITE 

FORMAT 


120 
90 3 
130 

140 


+ 

+ 

+ 

DO 120 
WRITE 
FORMAT 
DO 14 0 
KOUNL ( 
LOO PN ( 
KK = -K 
GO TO 


EQ . 1 
(4,9 
( 1 H 0 

1H 


STEP 6. CHECK FOR SINGULAR JACOBIAN 


ISIMDD, 


1 = 1 

(4,9 
( 1H 
1 = 1 
r ) =o 
i ) =o 

210 


******************* 

) GO TO 150 

02) K, ID, KOUNL( ID ) 

, ' ERROR IN NLESTN , MATRIX FAILURE USING 
' ERROR INDICATOR IS',13,/, 

' ' FIRST ARGUMENT IN THE CALL STATEMENT IS' 15 
' LOOP COUNTER IS ',13, 

THE COLUMN AND SQUARE MATRICES FOLLOW' ) 

, N 

03 ) B(I) , { A ( X , J ) , J»1 ,N) 

, G1 4 . 6 , 5X , 6 ( G1 4 . 6 ) ) 

, 10 


C 

c * * * 

150 


160 

170 

C 

(2 * * * 
18 0 

190 


7 - INCRE ^ENT INDEPENDENT VARIABLES WITHIN TRUST REGION 
DO 16 0 1=1 , N 

TF ACT= DABS ( D ( I , 1 ) ) / ( PCTMAX *DABS{SX(I) ) ) 

IF ( TFACT . GT . FACT ) F ACT=TF ACT 
DO 170 1=1 , N 
SX ( I ) =SX ( I } -D ( I , 1 ) /FACT 
L= 0 

GO TO 190 

8 ' R£SET ARGUMENT list independent variables and counters 

S X ( L ) = 5 X ( L ) * S F ( L ) 

KK= L+ 1 


X1=SX ( 1 ) 
X 2 = S X ( 2 ) 
X 3 = SX ( 3 ) 


X 4 — SX { 4 ) 

X 5 — SX ( 5 ) 

X 6 = S X ( 6 ) 

C 

200 LOO PN ( I D ) = L 
C 

210 RETURN 
END 
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Appendix 3. Subroutine NLESTN implementation with trust 
region Broyden Rank One method. 


♦ A******************, 


CODE DESIGNATION BROYDEN ********************** 
SUBROUTINE NLESTN {I D , N # A , KK , XI , Y1 , T1 , F 1 , X 2 , Y 2 T2 F2 X3 Y3 T3 FT 

> X4,Y4,T4,F4,X5,Y5,T5,F5,X6,Y6,T6,F6) ' ' ' ' # 

IMPLICIT RE AL * 8 { A-H , O - Z ) 


—NLESTN 2- TO 6 -DIMENS IONAL BROYDEN'S (GOOD) RANK ONE METHOD FOR 
THE ITERATIVE SOLUTION OF SIMULTANEOUS NONLINEAR ALGEBRAIC 
EQUATIONS . 


C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


CALL STATEMENT. . . 

CALL NLESTN (ID,N,A,KK,X1,Y1,T1,F1,X2,Y2,T2,F2. . . -XN f YN,TN FN ) 
. . .WHERE N IS >a6. ' ' ' 


ID IS THE ITERATION LOOP NUMBER (MINIMUM VALUE * 1, MAXIMUM 
VALUE =10). EACH NUMBER OF A NEST OF LOOPS MUST HAVE A DIFFERENT 
LOOP NUMBER, BUT OTHERWISE LOOP NUMBERS ARE ARBITRARY WITHIN THE 
ABOVE LIMITS. 

N IS THE NUMBER OF INDEPENDENT VARIABLES TO BE ITERATED UPON 
<X1,X2...XN) IN ORDER TO DIMINISH THE DEPENDENT VARAIBLt ERROR 
VALUES (Y1,Y2...YN) TO WITHIN SPECIFIED TOLERANCES (T1,T2 TN ) 
NOTE- 2 < = N < = 6 ' } # 

<F1,F2...FN) ARE INDEPENDENT VARIABLE INCREMENT MULTIPLIERS 
USED IN THE FINITE DIFFERENCE APPROXIMATION TO THE JACOBIAN 
MATRIX AT SELECTED STAGES. 

A rs AN N-BY-N MATRIX DIMENSIONED IN THE CALLING PROGRAM. 

KK IS A FLAG AS FOLLOWS... 

K K = - 3 SINGULAR MATRIX. 

KK=-t ALLOWABLE NUMBER OF ITERATIONS HAS BEEN EXCEEDED. 

KK= 0 ALL Y'S ARE LESS THAN TOLERANCE, ITERATION COMPLETE. 

KK= 1 ONE OR MORE Y'S ARE GREATER THAN TOLERANCE, REITERATE 


INCLUDE (INSAVE) 

INCLUDE ( DP SAVE ) 

DIMENSION NBROY (10) ,SX(6),SY0{10,6) ,SY(6) ,ST(6),SF(6) ,SDF(6) 
DIMENSION A( N, N ) , C ( 36 ) , SK ( 10 , 6 ) , SID ( 6 , 6 ) 

DATA LIMIT , PCTMAX , SMNUM/50 ,0.2,1. OD-8 / 

I.* *** V. ************* * ARRAY INITIALIZE SUBROUTINE INPUT ARGUMENTS 

S >. 1 ■ = N I 

s:-: f 2 , 
s :*: * 3 > = x 3 
5X ■ 4 : = X 4 
s>: ( 5 - =X 5 


SY ( 1 ■■ = 71 
SY( 2 }=Y2 
SY( 3 ) = Y 3 
SY( 4 ) = Y 4 
SV { 5 ) =Y5 
SY ( 6 )=Y6 


ORIGINAL PAGE IS 
OF POOR QUALITY 
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ST ( 1 ) =T 1 
ST ( 2 ) =T 2 
ST ( 3 )*T3 
ST ( 4 ) =T 4 
ST ( 5 ) = T 5 
ST( 6 )aT6 

ORIGIN^ 
OF POOH 

QUALITY 

SF ( 1 ) =*F 1 
SF( 2 }=F2 
SF(3 )=F3 
S F { 4 ) s=F 4 
SF( 5)=F5 
SF ( 6 ) s F 6 




10 


---...rw******* STE 

L = LOO PN (ID) 

NUMN ( ID, X )=NUMN( ID, 1 )+l 
I CONV^O 
DO 10 I a 1 , N 

IF (ICONV.GT.O) GO TO 10 

IF ( DABS ( S Y ( I ) ) . GT . DABS ( ST ( I ) ) ) ICONV=l 
CONTINUE i 

IF (ICONV.GT.O ) GO TO 20 

NUMN ( ID , 3 ) aNUMN ( ID , 3 } +1 

NUMN{ ID, 4 )=MAX0 ( NUMN ( ID, 4 ) , KOUNL < ID) I 

KOUNL(ID)=0 ' M 

KK = 0 

GO TO 200 


* * 

******** 

******* 

20 

IF { KOUNL { I D ) 


WRITE 

(4,901) 

9 0 1 

format 

( 


+ 1H , 

' ERROR 


+■ 1H p 

'FIRST . 


+■ 1H , 

'ERROR i 


+ 

K=1 

3 ( G i S . 6 


GO TO 

L 3 0 


^ LT ^LIMIT ^ *<30 TO E 30 2 * F ° R EXCESSIVE DERATIONS 
LIMIT, ID, ( SY ( I ) , ST ( I ) ,SX(I) , 1=1 , N ) 

:N NLESTN, no SOLUTION WITHIN ',13,' ITERATIONS ' / 
ARGUMENT IN THE CALL STATEMENT IS ',15 / ' A 
'ALUES TOLERANCES INDEPENDENT VARIABLES' ,/, 


c 

c 

c 

c 


30 


********, 


STEP 3 


IF ( KOUNL { ID } . LE . 0 ) 
IF ( NBROY ( ID) .GT . 0 ) 


NBROY ( ID) 
GO TO 300 


= 0 


PROCEDURE FOR 


DECIDE ON ITERATION 
CURRENT PHASE 
NBROY ( ID)=s0 NEWTON-RAPHSON STEP 
NBROY (ID) > 0 BROYDEN RANK 1 STEP 


4 0 


DO 4 0 

1 = 1 

SDF { I 

) = SF 

IF (L. 
L= 0 

, GT . 

DO 50 

1 = 1 


STEP 4. INITIALIZE ITERATION SEQUENCE FOR 
MULTIVARIATE newton-raphscn STAGE 


0 ) GO TO 60 
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OF POOR QuALfTv 


50 SYO ( ID, I )*SY ( I ) 

KOUNL { I D ) = KOUNL { I D ) + 1 
NUMN ( ID, 2 )=NUMN( ID, 2 )+l 
GO TO 180 


********************************** STEP 5 
60 SX(L)*SX(L)/SF(L) 

DP=SX(L> * S DF ( L) 

DO 100 1 = 1, N 

100 A( I , L ) * ( S Y ( I)— SYO (ID, I) )/DP 
IF (L.LT.N) GO TO 180 


ESTIMATE JACOBIAN PARTIALS 


*************************** STEP 6 


S CALE = 0 . 

K=ISIMDD ( N , N , -N , A , S I D , SCALE , C ) 


DETERMINE INVERSE JACOBIAN MATRIX 
MATRIX [A] AS RETURNED FROM 
FUNCTION ROUTINE ISIMDD IS THE 
INVERSE JACOBIAN APPROXIMATION 


* * * 


902 


120 

903 

130 

140 


+ 

+ 

+ 

+ 


* * * 


******** 


**************** STEP 7 


CHECK FOR SINGULAR JACOBIAN 


IF ( K.EQ . 1 ) GO TO 150 
WRITE (4,902) K , I D , KOUNL ( I D ) 

FORMAT ( 1H0 , ' ERROR IN NLESTN , MATRIX FAILURE USING ISIMDD,', 
' ERROR INDICATOR IS', 13,/, 

1H , 'FIRST ARGUMENT IN THE CALL STATEMENT IS ',15, 

' LOOP COUNTER IS ' , 13 , 

' THE COLUMN AND SQUARE MATRICES FOLLOW' > 


DO 120 1 = 1, N 

WRITE (4,903) SYO ( ID, I ) , ( A( I , J ) , J=1 , N ) 
FORMAT (1H ,G14.6,5X,6(G14.6) ) 

DO 140 1=1,10 
KOUNL ( I ) =0 
LOO PN ( I ) =0 


K K = - K 


GO TO 210 
C 

c ********* STEP 8. INCREMENT INDEPENDENT VARIABLES WITHIN TRUST REGION 
150 DO 155 1=1, N 
S K ( I D , I 1=0.0 
DO 155 J=1,N 

155 SK{ ID, I j = SK( ID, I )-A{ I , J) *SY0 ( ID,J ) 

F ACT = 1 . 0 
DO 160 1=1 ,N 

T FACT = DABS (SK(ID,I)/( PCTMAX*SX ( I ) ) ) 

160 IF ( TF ACT . GT . FACT ) F ACT = TF ACT 
DO 170 1=1 , N 

SK( ID, I ) =SK ( I D , I ) /FACT 
170 SX(I,'=SX'I)+SK(ID,I) 

L=0 


NBROY ( ID ) =1 
GO TO 190 
C 

C ************************* STEP 9 
C 

300 KOUNL ( ID )=KOUNL ( ID ) + l 


PROVIDE BROYDEN RANK O^E UPDATE FOR 
INVERSE JACOBIAN APPROXIMATION 
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ORIGINAL page is 

~ ~ j- ! / 

Of 


310 

320 


33 0 


340 

350 

355 


3 6 0 


370 


380 


NUMN{ ID, 2 )=NUMN( ID, 2 ) +1 
SBY=Q. 

DO 320 1=1 

S B = 0 . 


DO 310 J=1,N 
S B»S B + SK { ID, J) * A { J , I ) 

SBY-SBY + SB* <SY( I >-SY0 ( ID, I ) ) 

IF {DABS(SBY) . LT . SMNUM ) GO TO 355 
DO 350 1*1, N 
BY = 0 . 


DO 3 3 0^ 11*1 , N 

BY*BY+A( I,II)*(SY(II)-SY0(ID.II)) 
DO 350 J»1 ,N 
SB*0 . 

DO 340 JJ= 1 , N 
SB*SBvSK(ID,JJ)*A(JJ,J) 

A ( I , J ) * A (I,J) + (SK(ID,I)-BY) * SB/S BY 
DO 360 1*1, N 
SK( ID, I )*0 . 

DO 360 J*1,K 

SK( ID, I )-SK(ID f I J-A(I, JJ *SY( JJ 

F ACT* 1 . 0 
DO 370 1=1, N 

TF ACT* DABS < SK{ ID, I ) / ( PCTMAX * SX ( I ) ) ) 

IF { TFACT . GT . FACT ) F ACT=eTFACT 

DO 380 1=1, N 

SK( ID, I )*SK(ID,I ) /FACT 

SX ( I ) =SX ( I ) +SK ( ID , I ) 

SY0iID,I)=SY(I) 

NBROY ( ID)=NBROY( ID }>1 
IFq(NBROY(ID) . GT . N ) NBROY(ID)=0 


GO TO 190 


STEP 10. RESET ARGUMENT LIST INDEPENDENT VARIABLES AND COUNTERS 

SXIL)=SX(L)*SF<L) 

190 KK = L +■ 1 


xi = s;c{ l ) 

X 2 = S X f 2 ) 

X 3 * S X ( 3 ) 

X 4 = S X { 4 ) 

X5 = S>: f 5 ) 

X 6 = S X ( 6 } 

r 

2 00 LOOPN f I D ) =L 
C 

210 RETURN 
END 
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